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A numerical technique is proposed for an efficient numerical determination of the average phase 
factor of the fermionic determinant continued to imaginary values of the chemical potential. The 
method is tested in QCD with eight flavors of dynamical staggered fermions. A direct check of the 
validity of analytic continuation is made on small lattices and a study of the scaling with the lattice 
volume is performed. 
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I. INTRODUCTION 

Lattice QCD simulations in presence of a finite den- 
sity of baryonic matter are hindered by the well known 
sign problem. Consider for instance the QCD partition 
function 

Z{fi,[i) = J ■DUe- SG M(detM[U,(j]) 2 



VUe- SG[U] \detM[U,fj]\ 2 ( 



i2l) 



(1) 



describing two flavors of quarks (or eight flavors in the 
case of staggered flavors) which are given an equal chem- 
ical potential /1: the determinant of the fermionic ma- 
trix M is in general complex (9 7^ 0) for fi 7^ and 
Monte Carlo simulations are not feasible. Various pos- 
sibilities have been explored to circumvent the prob- 
lem, like reweighting techniques [l|, H, the use of 
an imaginary chemical potential cither for analytic con- 
tinuation [J, H H, 0, d, H 03 or for reconstructing 
the canonical partition function [ll|, [l2|, G3 , Taylor ex- 
pansion techniques fLU [l5| and non-relativistic expan- 
sions m in m . 

The same is not true in the case of a finite isospin 
density, i.e. when quarks are given opposite chemical 
potentials. Indeed, due to the property det M[U, — [/,] = 
det M[U,fi]*, the partition function 



ZQi,-n) = J VUe- SG[u] \detM[U,n}\ 2 



(2) 



has a positive measure. That is also known as phase 
quenched QCD. The average value of the phase factor 
of the fermionic determinant, (e l2e )( Mi _ M ), where the in- 
dex indicates the partition function the expectation value 
refers to, gives a direct measurement of the severeness of 
the sign problem. (e l2e ) ~ will signal the stage at which 
the complex nature of the determinant will imply a signif- 
icant difference between finite baryonic density and finite 
isospin density, as well a poor reliability of reweighting 
techniques (see Ref. [l9[ and references therein). 

It clearly follows from Eqs. ([T]) and © that the average 
phase factor is the expectation value of the ratio of two 



determinants; it can also be expressed as the ratio of two 
partition functions: 



i28\ 



I detM(/x) 



det M(-ju) 



g (p, m) 

Z(jJ,,-fJL) 



(3) 



Its direct numerical computation reveals a difficult nu- 
merical task as the lattice volume V increases, since it 
involves the numerical evaluation of fermionic determi- 
nants. 

It has been proposed recently [|(| HH to study the an- 
alytic continuation of the average phase factor to imagi- 
nary values of the chemical potential 



A28\ 



I III 



det M{i\i) \ Z{in,in) 
det M (-»m) / ( itl - ifl ) -¥) 
J VUe~ SG ^ det M[U, ifi] det M[U, i/j] 
J VUe- s oVJ] detM[J7,z//] det M[U,-ifi] 



(4) 



where Z(ifj,,ifj,) and Z(ifj,,—ifi) are the analytic contin- 
uation of the partition functions at finite baryonic and 
isospin chemical potentials respectively, which are both 
suitable for numerical simulations since det M[U,i/j] is 
always real. Numerical difficulties however are present 
also in this case: the observable to be averaged is still 
expressed in terms of fermionic determinants. Moreover 
in principle sampling problems deriving from a bad over- 
lap between the two statistical distributions described 
by Z(ifi,ifi) and Z(in,—ifi) may arise. In Ref. [2l| the 
fermionic determinant has been estimated on the basis of 
the lowest lying eigenvalues of the fermionic matrix. 

In the present paper we propose a new technique 
which, making use of numerical strategies developed in 
different contexts, permits an exact evaluation of the av- 
erage phase factor with a reasonable scaling of the re- 
quired CPU time as the lattice volume is increased. In 
doing this we will fully exploit the possibility of per- 
forming numerical simulations of the partition function 
Z{i\x\,i\i-i) for generic values of [i\ and [12- 

In Section [TT] we illustrate two different possible meth- 
ods, which are then numerically tested and compared in 
Section UTH for the theory with 8 staggered flavors. 
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II. THE METHOD 

The evaluation of the average phase factor, expressed 
like in Eq. or Eq. (TJJ as the ratio of two different par- 
tition functions, resembles similar problems which are 
encountered in quite different contexts, like the evalua- 
tion of disorder parameters in statistical models and in 
lattice gauge theories. Explicit examples are given by 
monopole disorder parameters or by the 't Hooft loop, 
which enters in various studies about color confinement. 
The major problem in those cases is the small overlap be- 
tween the statistical distributions corresponding to two 
different partition functions, resulting in a poor sampling 
efficiency. Powerful techniques have been developed in 
both cases, consisting in either determining derivatives 
of the disorder parameters, from which the ratio of par- 
tition functions can then be reconstructed after integra- 
tion [13, [H, HH, or in making use of various reweight- 
ing techniques, like that of rewriting the original ratio in 
terms of intermediate ratios which are more easily evalu- 

able [H IMIiS- 
In the present case the major difficulty derives from a 
direct computation of the observables, which is expressed 
in terms of fermionic determinants, but sampling prob- 
lems may in principle worsen the situation also in this 
case, especially in the large volume limit. In the fol- 
lowing we will describe the application of both kind of 
techniques described above to the present case, and try 
to understand by numerical simulations which of them is 
best suited in this context. 

We describe at first how to reconstruct (e l2e ) in terms 
of derivatives. Consider the modified ratio 



Z(ip, iv) 
Z(ifi, — ip,) 



(5) 



where 



Z(ifi,iv)= / VUe- Sc[u] det M[U,ip]det M[U,iv] . (6) 



It is clear that R^(—fi) = 1, while i? Al (/i) is the original 
ratio. It can be easily verified that 

pty) = ^lni? M (^) = -^-lnZ(i/i, ii/) 

•^'-V)^M(„,) )^. (7) 

Last quantity is nothing but i times the average num- 
ber of quarks coupled to the chemical potential iv: it 
is purely imaginary for symmetry reasons, hence p{y) is 
real, and can be computed using a noisy unbiased esti- 
mator. The average phase factor can then be obtained 
by integration 



(e l2 % = exp ^ p(v)d^j 



(8) 



In practice, the derivative p(y) will be computed for 
a discrete set of values of v and then integrated numeri- 
cally. The precision attained for (e l26 )i^ will depend both 
on the statistical errors of the single determinations and 
on the systematic uncertainty linked to numerical inte- 
gration; the last can be estimated for instance by varying 
the chosen interpolation procedure. In principle it is also 
possible to determine further derivatives of p in order to 
improve the integration accuracy. 

As a different method we consider rewriting (e 
the product of N intermediate ratios: 



■20\ 



as 



i28\ 



Z(ip, ip) 



N 



Z(ip, —ip) Zn 



Zn Zff—l Z\ -pi- 

.—=_[_[ r k (9) 



-1 



k=l 



where Zn = Z(ip,ip), Zq = Z(ip, —ip) while 



VUe- Salu] det M[U, ip] det M[U, i(-p + 1e5v)](lQ) 



Z, 



with Sv — 2p/N. The idea is to compute each single 
ratio rfc by a different Monte Carlo simulation. Apart 
from the increased overlap among each couple of partition 
functions, an improvement comes also from the simpler 
form in which the observable appearing in each ratio r k 
can be rewritten, for large enough N. Indeed we have: 

r fc = (det M(i(v + Sv))/ det M(iv)) {ifiiu) 

= (exp(TrlnA(v,6v))) {ifl . v) (11) 



where v = —p + (fc — l)5v and 

A[U, v, Sv] = M[U, iv]~ x M[U, i(y + 8v)\ . 



(12) 



If bv is small, the matrix A[U, v, 8v] is very close to 
the identity matrix Id for each configuration U. We can 
therefore expand the logarithm in Eq. pip thus rewriting 
the following approximate expression for r^: 



n , - { « xi. ( Tr(A - Id) - ~Tr(A - Id) 2 



(13) 



Each trace in the exponential can be evaluted by a 
noisy estimator as follows: 



1 K 

Tt(A[U] - Id)" ~ —J2 v {j)i ( A [u] - W l v 



(14) 



3=1 



and no quark determinant must be explicitly computed. 



where rf^ is a random vector satisfying (j^ ftf )n = 
Si u i 3 . The computation of each noise estimate in 
Eq. p4[) can be made faster if, when applying the ma- 
trix A[U] = M[U, iv]~ 1 M[U, i(v + 5v)] to the vector r]^ 
(or to (A[U] — Id) n ~ 1 r]^ at higher orders), rf^> itself is 
taken as a starting tentative solution for the inverter giv- 
ing M[U,wY 1 {M[U, i(y + Sv)]rj^): the guess is better 
and better as 5v — > 0. 

The second method is not conceptually different from 
the first one: finite free energy differences are computed 
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FIG. 1: p(v) for various values of p at /3 = 4.8 and L s = 4. 



FIG. 3: p(v) for various values of p at (3 — 4.8 and L a — 16. 



in this case instead of derivatives. However the numeri- 
cal procedure is different and it is not clear apriori which 
approach is more convenient. In the second case no nu- 
merical integration must be performed, however one has 
the drawback that the exponential of a noisy unbiased es- 
timator is biased, hence a large number K of vectors must 
be used and the final result must be checked to be inde- 
pendent of K. Moreover the systematic error involved 
in the truncation of the logarithm expansion, Eq. (p2 
must be properly estimated and kept under control. 



III. NUMERICAL RESULTS 

We have tested our methods for the theory with 8 stag- 
gered flavors of mass am = 0.1. We will present results 
obtained on x L t lattices with L t = 4 and L s = 4, 8, 16. 
At zero chemical potential this theory presents a strong 
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first order deconfinement/chiral transition, the critical 
coupling being (3 C ~ 4.7 for L t — 4. We have performed 
simulations both in the deconfined region ((3 — 4.8) and 
in the low temperature confined region (/3 = 4.6). On 
the smallest lattice (L s = 4) we will compare our results 
directly with those obtained at real isospin chemical po- 
tential by a direct evaluation of the determinant phase. 
Numerical simulations have been performed mostly on 
the APEmille facility in Pisa. The INFN apeNEXT fa- 
cility in Rome has been used for the results on the largest 
lattice. The standard exact HMC algorithm [28| has been 
used with trajectories of length 1. 

A collection of the results obtained for imaginary 
chemical potentials is reported in Table |TJ 



1.06 




FIG. 2: p{v) for various values of \jl at /3 = 4.6 and L s — 4. 



FIG. 4: r~l (blank triangles) and the inverse of r^T (blank cir- 
cles) defined in Eq. (|15[1 . together with their geometric mean, 
i.e. ^pFk (filled circles). Data are showed for L 3 = 4, p = 0.05, 
v = -0.045 and 8u = 0.005. 
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A. Systematic errors and comparison of the 
methods 

In Figs. [T] [2] and [3] we report various determinations 
of p(y) (minus the imaginary part of the baryon number, 
Eq. {?])) obtained on discrete sets of points. 

It is apparent that p{y) is a very smooth function of 
v in all explored cases and independently of the lattice 
size and of the explored phase (confined or deconfined) 1 . 
In most cases it can even be approximated by a linear 
function; therefore numerical integration turns out to be 
an easy task. We have adopted a simple linear interpo- 
lation between consecutive points to obtain the results 
given in Table [H the reported errors derive from stan- 
dard error propagation of the statistical errors of the sin- 
gle data points. We have verified, by changing the order 
of the interpolating polynomial, that the systematic er- 
ror related to the interpolation-integration procedure is 
negligible with respect to the statistical one. 

Concerning the second method described in Section II, 
we have adopted a standard trick (2f| in order to reduce 
systematic effects. Each partial ratio rfc in Eq. (J9j> has 
been rewritten as 

_ / det M(i(v + 6v)) \ _r£_ 

_ (det Mtyy + 5v))l det Mfty + 

(detM(zH)/detM^+f))) v ^ +4f) (15) 

rfc can again be evaluated in a single simulation and a 
jackknife analysis has to be applied to obtain a correct 
error estimate. Two major benefits derive in this case. 
First, the reduced value 8v/2 greatly improves the con- 
vergence of the logarithm expansion in Eq. (|13|) . Second, 
the bias introduced by the finite number of noisy esti- 
mators, see Eqs. (Ti~3")) and (TT4")) . gets largely cancelled in 
the ratio. That is apparent from Fig. [4J where we plot 
ft and the inverse of defined in Eq. (|15p . and their 
geometric mean (i.e. y/Fk), as a function of the number 
K of noise vectors, in one particular sample case. It is 
clear that, while the single factors have a relatively slow 
convergence, their product is stable from K = 5 on. We 
have however always used K = 30 in our determinations. 

Regarding the logarithm expansion, Eq. (fT3|) . we have 
always adopted a third order approximation: in all cases 
the discrepancy with the result obtained at the second 
order is at least one order of magnitude smaller than the 
statistical uncertainty. The fact that the systematic error 
related to this expansion is well under control can be also 
appreciated from Table Q] second and third row, showing 



Clearly one expects a non-smooth behaviour if Z(i/x, i/x) and 
Zq = Z(ifi, —ifi) belong to two different phases, so that some 
transition is met when v goes from —v to v, however in these 
cases analytic continuation itself is not applicable. 
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TABLE I: Collection of determinations of the average phase 
factor continued to imaginary values of /i for various param- 
eter sets and computation methods. In the fourth column 
the method used to obtain the determination is described: 
DER(N) stands for the integration of the first derivative p de- 
termined on a discrete set of (N+l) points, Eq. JHJ; RAT(N) 
stands for the evaluation of N intermediate ratios rt , Eq. ((9} . 
Finally on the smallest lattices also a direct determination of 
the expectation value in Eq. (|4]) is reported for comparison. 



that the determination of (e l2e )i M is stable against the 
variation of the number of intermediate ratios. 

Let us now come to the comparison between the two 
methods. While they always give perfectly compatible 
results, thus confirming the absence of appreciable sys- 
tematics, it is clear from Table Q] that, with a compa- 
rable numerical effort (in the last column we give the 
total number of Monte-Carlo trajectories used for each 
determination), the method described by Eq. ([5]) (inte- 
gration of the derivative) furnishes more accurate deter- 
minations. We have therefore chosen this method in or- 
der to perform more extensive studies of (e 



B. Test of analytic continuation 

The average phase factor computed at finite isospin 
chemical potential, at variance with that computed in 
the quenched theory, is expected 0, Hl| to be an ana- 
lytic function of /i 2 around /i 2 = 2 . We can test directly 
analytic continuation by comparing our results with di- 
rect determinations of (e l2e )^ performed at real chemi- 
cal potentials: this is done only for the smallest lattice 
[L s = 4), where the second determination is easily af- 
fordable. 



2 It is even in fi for symmetry reasons. 
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FIG. 5: (e !2S ) computed for different values of /j 2 at /3 = 4.8 
and /? = 4.6 on a 4 4 lattice. Best fit quadratic and quartic 
functions in /x 2 are displayed, showing good validity of ana- 
lytic continuation. 



as shown in the figure. We obtain, at j3 = 4.8, A = 
-4.48(8), B = 15.7 ± 2.5 and x 2 /d.oi. ~ 1.3. 

Analyticity around /i 2 = is therefore well verified. 
We stress that at f3 = 4.8 our largest value of the imag- 
inary chemical potential is still below the first Roberge- 
Weiss phase transition at Im(/i) = 7r/(3L t ), hence within 
the expected range of validity of analytic continuation for 
/j 2 < at high temperatures. 



C. Large volume scaling 

We have performed numerical simulations at different 
values of L s in order to test both the behaviour of (e 120 ) 
and the efficiency of our method as the lattice volume is 
increased. 

In Fig. [S] we report determinations performed at fixed 
values of ifj, and variable L s at (3 = 4.8. A behaviour 

(O = 1 + CL1 (18) 



We plot in Fig. [5] results obtained at j3 = 4.6 and 
j3 = 4.8. The whole set of results obtained at real chem- 
ical potentials (/i 2 > 0) and imaginary chemical poten- 
tials (/i 2 < 0) can be described by a simple quadratic 
behaviour 

(e l2e ) = 1 + Afi 2 (16) 

in a range |/i 2 | < 0.01, with A = -4.41(9) and 
X 2 /d.o.f. ~ 1.5 for f3 = 4.8 and A = -10.2(3) and 
X 2 /d.o.f. ~ 1.8 for (3 — 4.6. If the range of values is 
extended a quartic term is necessary 

(e l2e ) - 1 + V + B M 4 (17) 



well describes the data with 7 ~ 2.5 for both values of 
'/'• 

Concerning the numerical efficiency, we notice that to 
obtain comparable uncertainties (of the order of 10 %) 
for (e l2e ) — 1, on the largest lattice (16 3 x 4) we needed a 
CPU time which is less than one order of magnitude big- 
ger than what needed on the smallest lattice (4 4 ). Con- 
sidering that the two lattice volumes differ by a factor 64, 
we deduce that, at least for the quark mass considered in 
the present study, our method requires a numerical effort 
which scales in an affordable way with the lattice size. 




FIG. 6: (e i2e ) ip. as a function of the spatial lattice size L s for 



two values of A best fit according to Eq. 
in both cases. 



18)1 is reported 



IV. CONCLUSIONS 

We have presented two different techniques, described 
respectively by Eq. (|5|) and Eq. ©, for an efficient nu- 
merical determination of the average phase factor of the 
fcrmionic determinant continued to imaginary values of 
the chemical potential. We have applied both methods 
to QCD with 8 dynamical staggered flavors, verifying the 
absence of uncontrolled systematic effects and perform- 
ing a comparison of the efficiencies, with the conclusion 
that the method based on the integration of the imagi- 
nary part of the baryon density, Eq. ([5]) , is numerically 
more convenient. A fair good scaling of the efficiency is 
observed as the lattice volume is increased. We have also 
directly tested, on small lattices, the analiticity of the av- 
erage phase factor around = 0. The method proposed 
and tested in the present paper will be used in the future 
to perform more extensive studies, with more physical 
quark masses and number of flavors, of the average phase 
factor continued to imaginary chemical potentials. 
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